Linear elasticity

Wrap-up: notebook

Author

Affiliation

Julia Kowalski

Chair of Methods for Model-Based Development in Computational Engineering

1D wave equation

\[\frac{\partial^2}{\partial t^2} \mathbf w = c \, \frac{\partial^2}{\partial x^2} \mathbf w\]

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

%matplotlib inline

1. Numerical solver based on finite differencing strategy:

def solve_wave_equation(c, L, T, dx, dt, initial_condition):
    
    # check stability
    if c * dt / dx > 1:
        raise ValueError("Stability condition violated: c*dt/dx should be <= 1")
    
    # discretize space-time domain
    x = np.arange(0, L + dx, dx)
    t = np.arange(0, T + dt, dt)
    u = np.zeros((len(t), len(x)))
    
    # set initial condition
    u[0, :] = initial_condition(x)
    if len(t) > 1:
        u[1, :] = u[0, :]  
    
    # buffer coefficient 
    C2 = (c * dt / dx)**2
    
    # time update
    for n in range(1, len(t) - 1):
        for i in range(1, len(x) - 1):
            u[n + 1, i] = (2 * u[n, i] - u[n - 1, i] + C2 * (u[n, i + 1] - 2 * u[n, i] + u[n, i - 1]))
    
    # set boundary condition
    u[:, 0] = 0
    u[:, -1] = 0
    
    return x, t, u

2. Set up simulation scenario

# Initial condition: Gaussian bump
def initial_gaussian_pulse(x, sigma=0.1, mu=0.5):
    return np.exp(-((x - mu) / sigma)**2)

# Scenario parameters
sigma = 0.25 # IC dev
mu = 0.5 # IC center
c = 1.0  # Wave speed
T = 2.0  # Total time
L = 1.0  # Domain size
dx = 0.01  
dt = 0.005

3. Deploy scenario 1 and 2

# Run simulation
x, t, u = solve_wave_equation(c, L, T, dx, dt, lambda x: initial_gaussian_pulse(x, sigma, mu))

# Visualize result
plt.figure(figsize = (10,8))
plt.imshow(u, extent=(0, L, 0, T), aspect='auto', origin='lower')
plt.colorbar(label='Wave amplitude')
plt.ylabel('Time')
plt.xlabel('Position')
plt.title('Wave Propagation')
plt.show()

Figure 1: Space-time plot at wave speed 1

# Run simulation
c = 0.5
x, t, u = solve_wave_equation(c, L, T, dx, dt, lambda x: initial_gaussian_pulse(x, sigma, mu))

# Visualize result
plt.figure(figsize = (10,8))
plt.imshow(u, extent=(0, L, 0, T), aspect='auto', origin='lower')
plt.colorbar(label='Wave amplitude')
plt.ylabel('Time')
plt.xlabel('Position')
plt.title('Wave Propagation')
plt.show()

Figure 2: Space-time plot at wave speed 0.5

4. Animate result

%%capture
# Run simulation
c = 1.0
x, t, u = solve_wave_equation(c, L, T, dx, dt, lambda x: initial_gaussian_pulse(x, sigma, mu))


def animate_wave(x, u, interval=50):
    fig, ax = plt.subplots()
    line, = ax.plot(x, u[0, :], 'k-')
    ax.set_xlim(x.min(), x.max())
    ax.set_ylim(-1.1, 1.1)
    ax.set_xlabel('Position')
    ax.set_ylabel('Wave amplitude')
    ax.set_title('1D Wave Propagation')

    def update(frame):
        line.set_ydata(u[frame, :])
        return line,

    ani = FuncAnimation(fig, update, frames=len(u), interval=interval, blit=True)
    #plt.close(fig)
    return ani

ani = animate_wave(x, u)

HTML(ani.to_html5_video())
# ani.save('./WavePropagation.mp4', writer='ffmpeg', dpi=300)

# Run simulation
c1 = 1.0
c2 = 0.5
x, t, u1 = solve_wave_equation(c1, L, T, dx, dt, lambda x: initial_gaussian_pulse(x, sigma, mu))
x, t, u2 = solve_wave_equation(c2, L, T, dx, dt, lambda x: initial_gaussian_pulse(x, sigma, mu))

# Visualize result
plt.figure(figsize = (10,8))
# fig, ax = plt.subplots(figsize=(10,8))
plt.imshow(u1+u2, extent=(0, L, 0, T), aspect='auto', origin='lower')
plt.colorbar(label='Wave amplitude')
plt.ylabel('Time')
plt.xlabel('Position')
plt.title('Wave Propagation')
plt.show()

Figure 3: Space-time plot for two superposed waves

%%capture

u_superposed = 0.5*(u1 + u2)

def animate_wave(x, u, interval=50):
    fig, ax = plt.subplots()
    line, = ax.plot(x, u[0, :], 'b-')
    ax.set_xlim(x.min(), x.max())
    ax.set_ylim(-1.1, 1.1)
    ax.set_xlabel('Position')
    ax.set_ylabel('Wave amplitude')
    ax.set_title('1D Wave Propagation')

    def update(frame):
        line.set_ydata(u[frame, :])
        return line,

    ani = FuncAnimation(fig, update, frames=len(u), interval=interval, blit=True)
    #plt.close(fig)
    return ani

ani = animate_wave(x, u_superposed)

HTML(ani.to_html5_video())
# ani.save('./WavePropagationSuperposed.mp4', writer='ffmpeg', dpi=300)

%%capture


def animate_wave(x, u, u_superposed, interval=50):
    fig, ax = plt.subplots()
    line, = ax.plot(x, u[0, :], 'k-', label='single wave')
    line2, = ax.plot(x, u_superposed[0, :], 'b-', label='superposed wave')
    ax.set_xlim(x.min(), x.max())
    ax.set_ylim(-1.1, 1.1)
    ax.set_xlabel('Position')
    ax.set_ylabel('Wave amplitude')
    ax.set_title('1D Wave Propagation')
    plt.legend()

    def update(frame):
        line.set_ydata(u[frame, :])
        line2.set_ydata(u_superposed[frame, :])
        return line, line2

    ani = FuncAnimation(fig, update, frames=len(u), interval=interval, blit=True)
    #plt.close(fig)
    return ani

ani = animate_wave(x, u, u_superposed)

HTML(ani.to_html5_video())
# ani.save('./WavePropagationBoth.mp4', writer='ffmpeg', dpi=300)